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on : The performance of Gallager's error-correcting code is investigated via methods of statistical 
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physics. In this approach, the transmitted codeword comprises products of the original message 
bits selected by two randomly-constructed sparse matrices; the number of non-zero row/column 
elements in these matrices constitutes a family of codes. We show that Shannon's channel capacity 
is saturated for many of the codes while slightly lower performance is obtained for others which may 
be of higher practical relevance. Decoding aspects are considered by employing the TAP approach 

■ which is identical to the commonly used belief-propagation-based decoding. 

VD 

The ever increasing information transmission in the modern world is based on communicating messages reliably 
through noisy transmission channels; these can be telephone lines, deep space, magnetic storing media etc. Error- 
^-j correcting codes play an important role in correcting errors incurred during transmission; this is carried out by 
encoding the message prior to transmission, and decoding the corrupted received codeword for retrieving the original 
message. In his ground breaking papers, Shannon jjj analyzed the capacity of communication channels, setting an 
upper bound to the achievable noise-correction capability of codes, given their code (or symbol) rate. The latter 
represents the ratio between the number of bits in the original message and the transmitted codeword. 

Shannon's bound is non-constructive and does not provide explicit rules for devising optimal codes. The quest for 
more efficient codes, in the hope of saturating the bound set by Shannon, has been going on ever since, providing 
many useful but sub-optimal codes. 

• One family of codes, presented originally by Gallager attracted significant interest recently as it has been shown 
to outperform most currently used techniques In fact, irregular versions of Gallager-type codes have recently been 
shown to get very close to saturating Shannon's bound in the case of infinitely long messages [Q . Gallager-type codes 
are characterized by several parameters, the choice of which defines a particular member of this family of codes. Most 

■^j- ■ studies of Gallager-type codes conducted so far have been carried out via numerical simulations. Some analytical 
results have been obtained via methods of information theory |s| , setting bounds on the performance of certain code 

• types, and by combinatorical/statistical methods B; no quantitative results have been obtained for their typical 
| performance. 

In this Letter we analyze the typical performance of Gallager-type codes for several parameter choices via methods 
Q\ [ of statistical mechanics. We then validate the analytical solution by comparing the results to those obtained by the 

■ TAP approach to diluted systems and via numerical methods. 

In a general scenario, a message represented by an N dimensional Boolean/binary vector £ is encoded to the M 
dimensional vector J° which is then transmitted through a noisy channel with some flipping probability p per bit 
(other noise types may also be considered but will not be examined here). The received message J is then decoded 
to retrieve the original message. 

One can identify several slightly different versions of Gallager-type codes. The one used in this Letter, termed the 
O ■ MN code Q is based on choosing two randomly-selected sparse matrices A and B of dimensionality MxN and MxM 
respectively; these are characterized by K and L non-zero unit elements per row and C and L per column respectively. 
The finite, usually small, numbers K, C and L define a particular code; both matrices are known to both sender and 
receiver. Encoding is carried out by constructing the modulo 2 inverse of B and the matrix B~ 1 A (modulo 2); the 
vector J° =B~ 1 A £ (modulo 2, £ in a Boolean representation) constitutes the codeword. Decoding is carried out by 
taking the product of the matrix B and the received message J = J°+C (modulo 2), corrupted by the Boolean noise 
vector £, resulting in A£ + B£. The equation 

A£ + BC = AS + Bt (I) 

is solved via the iterative methods of Belief Propagation (BP) |1 to obtain the most probable Boolean vectors S and 
r; BP methods in the context of error-correcting codes have recently been shown to be identical to a TAP || based 
solution of a similar physical system 

The similarity between error-correcting codes of this type and Ising spin systems was first pointed out by Sourlas 
0, who formulated the mapping of a simpler code, somewhat similar to the one presented here, onto an Ising spin 
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system Hamiltonian. We recently extended the work of Sourlas, that focused on extensively connected systems, to 
the finite connectivity case 0. 

To facilitate the current investigation we first map the problem to that of an Ising model with finite connectivity. 
We employ the binary representation (±1) of the dynamical variables S and r and of the vectors J and J° rather 
than the Boolean (0, f ) one; the vector J° is generated by taking products of the relevant binary message bits 
J?.^ j;i . = ^i ± ^i 2 . . ., where the indices ■ ■ ■ correspond to the non-zero elements of B~ 1 A 1 producing a binary 
version of J°. As we use statistical mechanics techniques, we consider the message and codeword dimensionality 
(N and M respectively) to be infinite, keeping the ratio between them R = N/M, which constitutes the code rate, 
finite. Using the thermodynamic limit is quite natural as Gallager-type codes are usually used for transmitting long 
(10 4 — 10 5 ) messages, where finite size corrections are likely to be negligible. To explore the system's capabilities we 
examine the Hamiltonian 



H 



^ t T^<il,..,iK\jl,--,jL> & 



<i 1 ,..,i K ;j 1 ,..,j L > 

Ojj . . . Si K Tj 1 . . . Tj L 



N 



1 i J<ii,..,iK;ji,--jL> 
M 



(2) 



The tensor product V <iu , .,i Ki j 1 ,..j L >J<.i 1 ,..,i K ;j 1 ,..ji,>, where J<i u ..,j L> = • • • ZiicGjith • ■ ■ Ql > is the binary 

equivalent of A£+B£, treating both signal (S and index i) and noise (t and index j) simultaneously. Elements of the 
sparse connectivity tensor "D<i u ..,j L > take the value f if the corresponding indices of both signal and noise are chosen 
(i.e., if all corresponding indices of the matrices A and B are 1 ) and otherwise; it has C unit elements per i-index 
and L per j-index representing the system's degree of connectivity. The 8 function provides 1 if the selected sites' 
product Si 1 . . . Si K Tj 1 . . . Tj L is in disagreement with the corresponding element J<i x ,..j h >, recording an error, and 
otherwise. Notice that this term is not frustrated, as there are M+N degrees of freedom and only M constraints from 
Eq.(|l]), and can therefore vanish at sufficiently low temperatures. The last two terms on the right represent our prior 
knowledge in the case of sparse or biased messages F s and of the noise level F T and require assigning certain values to 
these additive fields. The choice of /3^oo imposes the restriction of Eq. ([[]), limiting the solutions to those for which 
the first term of Eq.(@) vanishes, while the last two terms, scaled with j3, survive. Note that the noise dynamical 

variables r are irrelevant to measuring the retrieval success m = -k (J^iLi £i sign(jSj) fl ) . The latter monitors 
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the normalized mean overlap between the Bayes-optimal retrieved message, shown to correspond to the alignment of 
{Si)p to the nearest binary value Q], and the original message; the subscript (5 denotes thermal averaging. 

Since the first part of Eq.(||) is invariant under the transformations Si — » S^, Tj — * t,-£j and J<i lt ..j L > — » 
J<ii,--,jL>£.ii-4iK(jiCj2--(jL — 1; it would be useful to decouple the correlation between the vectors S, r and £, ^. 
Rewriting Eq.(Q) one obtains a similar expression apart from the last terms on the right which become F s / P^2 k Sk £fc 
and F T /(3J2k T k Ck- 

The random selection of elements in T> introduces disorder to the system which is treated via methods of statistical 
physics. More specifically, we calculate the partition function Z(T>, J) = Tr^g x j exp[— (3H] averaged over the disorder 
and the statistical properties of the message and noise, using the replica method f]]|||. Taking /3^oo gives rise to 
a set of order parameters 
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where a, 0, .. represent replica indices, and the variables Zi and Yj come from enforcing the restriction of C and L 
connections per index respectively H: 
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and similarly for the restriction on the j indices. 
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To proceed with the calculation one has to make an assumption about the order parameters symmetry. The 
assumption made here, and validated later on, is that of replica symmetry in the following representation of the order 
parameters and the related conjugate variables 



la, p.. i — a q dx tt(x) x l , q a ,p..-y = I dx tt(x) x 
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r a ,f3..i = a r / dy p(y) y , f a ^.. 7 = a~ / dy p(y) y 



where I is the number of replica indices, a* are normalization coefficients, and 7r(x), tt (x), p(y) and p{y) represent 
probability distributions. Unspecified integrals are over the range [— 1,+1]. One then obtains an expression for the 
free energy per spin expressed in terms of these probability distributions 
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where (-}^ and (-)^ denote averages over the input and noise distributions of the form 



El 1 + tanh F s 1 - tanh F s . 
1 5 + o f (0 
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and similarly for (•)<; where is replaced by F T . 

The free energy can then be calculated via the saddle point method. Solving the equations obtained by varying 
Eq.(^|) w.r.t the probability distributions 7r(x), t?(x), p(y) and p(y), is generally difficult. The solutions obtained in 
the case of unbiased messages (the most interesting case as most messages are compressed prior to transmission) are 
for the ferromagnetic phase: 



n(x) — 5(x — 1) , 7? (x) = S(x — 1) 
p(y) = S(y - 1) , p(y) = S(y - 1) , 

and for the paramagnetic phase (there is no spin-glass solution due to lack of frustration): 



(8) 
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It is easy to verify that these solutions obey the saddle point equations. However, it is necessary to validate the 
stability of the solutions and the replica symmetric ansatz itself. To address these questions we obtained solutions 
to the system described by the Hamiltonian (||) via the TAP method of finitely connected systems H ; we solved the 
saddle point equations derived from Eq.(^j) numerically, representing all probability distributions by up to 10 4 bin 
models and by carrying out the integrations via Monte-Carlo methods; finally, to show the consistency between theory 
and practice we carried out large scale simulations for several cases, which will be presented elsewhere. The results 
obtained by the various methods are in complete agreement. 

The various methods indicate that the solutions may be divided to two different categories characterized by K = L = 2 
and by either K >3 or L>3, which we therefore treat separately. 

For unbiased messages and either K > 3 or L > 3 we obtain the solutions @ and ([)]) both by applying the TAP 
approach and by solving the saddle point equations numerically. The former was carried out at the value of F T which 
corresponds to the true noise and input bias levels (for unbiased messages F s = 0) and thus to Nishimori's condition 
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, where no replica symmetry breaking effect is expected. This is equivalent to having the correct prior within the 
Bayesian framework and enables one to obtain analytic expressions for some observables as long as some gauge 
requirements are obeyed pj. Numerical solutions show the emergence of stable dominant delta peaks, consistent with 
those of (||) and @ . The question of longitudinal mode stability (corresponding to the replica symmetric solution) 
was addressed by setting initial conditions for the numerical solutions close to the solutions ^ and (|^), showing that 
they converge back to these solutions which are therefore stable. 

The most interesting quantity to examine is the maximal code rate, for a given corruption process, for which 
messages can be perfectly retrieved. This is defined in the case of K, L>3 by the value of R = K/C = N/M for which 
the free energy of the ferromagnetic solution becomes smaller than that of the paramagnetic solution, constituting a 
first order phase transition. A schematic description of the solutions obtained is shown in the inset of Fig. la. The 
paramagnetic solution (m = 0) has a lower free energy than the ferromagnetic one (low/high free energies are denoted 
by the thick and thin lines respectively, there are no axis lines at m = 0, 1) for noise levels p>p c and vice versa for 
P<Pc] both solutions are stable. The critical code rate is derived by equating the ferromagnetic and paramagnetic 
free energies to obtain 

i2c = l-ffa(p) = l+(plogaP+(l-p)log a (l-p)) • (10) 

This coincides with Shannon's capacity. To validate these results we obtained TAP solutions for the unbiased message 
case (K = L = 3, C = 6). Averages over 10 solutions obtained for different initial conditions in the vicinity of the stable 
solutions are presented in Fig. la (as +) in comparison to Shannon's capacity (solid line). 

Analytical solutions for the saddle point equations cannot be obtained for the case of biased patterns and we 
therefore resort to numerical methods and the TAP approach. The maximal information rate (i.e., code rate x_ff 2 (/ s = 
(1 + tanhi 7 ' s )/2) - the source redundancy) obtained by the TAP method (O) and numerical solutions of the saddle 
point equations (□), averaged for each noise level over solutions obtained for 10 different starting points in the vicinity 
of the analytical solution, are shown in Fig. la. Numerical results have been obtained using 10 3 - 10 4 bin models for 
each probability distribution and had been run for 10 5 steps per noise level point. The various results are highly 
consistent and practically saturate Shannon's bound for the same noise level. 

The MN code for K,L > 3 seems to offer optimal performance. However, the main drawback is rooted in the 
co-existence of the stable m = 1 and m = solutions, shown in Fig. la (inset), which implies that from some initial 
conditions the system will converge to the undesired paramagnetic solution. Moreover, studying the ferromagnetic 
solution numerically shows a highly limited basin of attraction, which becomes smaller as K and L increase, while the 
paramagnetic solution at to = always enjoys a wide basin of attraction. As initial conditions for the decoding process 
are typically of close-to-zero magnetization (almost no prior information about the original message is assumed) it 
is highly likely that the decoding process will converge to the paramagnetic solution. This performance has been 
observed via computer simulations by us and by others || . 

While all codes with K, L > 3 saturate Shannon's bound and are characterized by a first order, paramagnetic to 
ferromagnetic, phase transition, codes with K = L = 2 show lower performance and different physical characteristics. 
The analytical solutions (j|) and @ are unstable at some flip rate levels and one resorts to solving the saddle point 
equations numerically and to TAP based solutions. The picture that emerges is sketched in the inset of Fig. lb: 
The paramagnetic solution dominates the high flip rate regime (appearing as a dominant delta peak in the numerical 
solutions) up to the point p\ (denoted as 1 in the inset) in which a stable, ferromagnetic solution, of higher free energy, 
appears (thin lines at to = ±1). At a lower flip rate valuer the paramagnetic solution becomes unstable (dashed line) 
and is replaced by two stable sub-optimal ferromagnetic (broken symmetry) solutions which appear as a couple of 
peaks in the various probability distributions; typically, these have a lower free energy than the ferromagnetic solution 
until pz, after which the ferromagnetic solution becomes dominant (at some code rate values it is dominant directly 
following the disappearance of the paramagnetic solution). Still, only once the sub-optimal ferromagnetic solutions 
disappear, at the spinodal point p 8 , a unique ferromagnetic solution emerges as a single delta peak in the numerical 
results (plus a mirror solution). The point in which the sub-optimal ferromagnetic solutions disappear constitutes the 
maximal practical flip rate for the current code rate and was defined numerically (O) and via TAP solutions (+) as 
shown in Fig. lb. 

Notice that initial conditions for both TAP and the numerical solutions were chosen almost randomly, with a 
very slight bias of 0(1O~ 12 ), in the initial magnetization. The TAP dynamical equations are identical to those 
used for practical BP decoding ||, and therefore provide equivalent results to computer simulations with the same 
parameterization, supporting the analytical results. The excellent convergence results obtained point out the existence 
of a unique pair of global solutions to which the system converges (below p s ) from practically all initial conditions. 
This observation and the practical implications of using the K = L = 2 code have not been obtained by information 
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theory methods (e.g. ||); these prove the existence of very good codes for C, L>3, and examine decoding properties 
only via numerical simulations. 

In this Letter we examined the typical performance of Gallager-type codes. We discovered that for a certain choice 
of parameters, either K > 3 or L > 3, one obtains optimal performance, saturating Shannon's bound. This comes at 
the expense of a decreasing basin of attraction making the decoding process increasingly impractical. Another code, 
K = L = 2, shows close to optimal performance with a very large basin of attraction, making it highly attractive 
for practical purposes. Studying the typical performance of Gallager-type codes, which complements the methods 
used in the information theory literature, is the first step towards understanding their exceptional performance and 
in the search for a principled method for designing optimal Gallager-type codes. Important aspects that are yet to 
be investigated include other noise types, irregular constructions and the significance of finite size effects. 
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FIG. 1. Critical code rate as a function of the flip rate p, obtained from numerical solutions and the TAP approach (N — 10 4 ), 
and averaged over 10 different initial conditions with error bars much smaller than the symbols size, (a) Numerical solutions 
for K = L = 3, C' — 6 and varying input bias f s (□) and TAP solutions for both unbiased (+) and biased (O) messages; initial 
conditions were chosen close to the analytical ones. The critical rate is multiplied by the source information content to obtain the 
maximal information transmission rate, which clearly does not go beyond 7? = 3/6 in the case of biased messages; for unbiased 
patterns H^ifs) — !- Inset: The ferromagnetic and paramagnetic solutions as functions of p; thick and thin lines denote stable 
solutions of lower and higher free energies respectively, (b) For the unbiased case of K = L = 2; initial conditions for the TAP 
(+) and the numerical solutions (O) are of almost zero magnetization. Inset: The ferromagnetic (optimal/sub-optimal) and 
paramagnetic solutions as functions of p; thick and thin lines are as in (a), dashed lines correspond to unstable solutions. 
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